!------------------------------------------------------!
!    Advection part   by Zhang Xuke     2021.9,7       !
!------------------------------------------------------!
        SUBROUTINE ADVECTION()
        USE BASIC_MOD
        IMPLICIT NONE
        INTEGER:: I,J,N,NN
        REAL*8:: FI(4)
        DO I=0,XCELLS
            DO J=0,YCELLS
                CALL AUSM(1,I,J,FI)
                FLUXE(I,J,:)=FI(:)
                FLUXW(I+1,J,:)=FI(:)
                DO NN=1,4
                IF(ISNAN(FI(NN))) THEN
                    WRITE(*,*) 'S',I,J
                    STOP
                END IF
                END DO
                CALL AUSM(2,I,J,FI)
                FLUXN(I,J,:)=FI(:)
                FLUXS(I,J+1,:)=FI(:)
                
            END DO
        END DO
        RETURN
        END SUBROUTINE ADVECTION
!----------------------------------------------------------------------! 
! PURPOSE:                                                             !
!----------------------------------------------------------------------!
      SUBROUTINE AUSM(DIRECTION,I,J,FI)
      USE BASIC_MOD
      IMPLICIT NONE
      INTEGER:: DIRECTION,I,J,SP
      REAL*8:: MAI,MAL,MAR,CI,CL,CR,VI,B1,G1,B2,G2,SIN,COS,VNL,VNR,PL,PR,DL,DR,VXL,VXR,VYL,VYR,HL,HR,GAMML,GAMMR,FI(4)
      IF(DIRECTION==1) THEN
	  IF(MUSCL.EQ.1) THEN
          CALL INTERPOLATION(PP(I-1:I+2,J),PL,PR,LIMIT) 
          CALL INTERPOLATION(DD(I-1:I+2,J),DL,DR,LIMIT) 
          CALL INTERPOLATION(UU(I-1:I+2,J),VXL,VXR,LIMIT) 
          CALL INTERPOLATION(VV(I-1:I+2,J),VYL,VYR,LIMIT) 
          CALL INTERPOLATION(HH(I-1:I+2,J),HL,HR,LIMIT) 
          CALL INTERPOLATION(GA(I-1:I+2,J),GAMML,GAMMR,LIMIT) 
	  ELSE 
          PL=PP(I,J)
	    PR=PP(I+1,J)
          DL=DD(I,J)
	    DR=DD(I+1,J)
          VXL=UU(I,J)
	    VXR=UU(I+1,J)
          VYL=VV(I,J)
	    VYR=VV(I+1,J)
          HL=HH(I,J)
	    HR=HH(I+1,J)
          GAMML=GA(I,J)
	    GAMMR=GA(I+1,J)
	  END IF 
      ELSE IF(DIRECTION==2) THEN
	  IF(MUSCL.EQ.1) THEN
          CALL INTERPOLATION(PP(I,J-1:J+2),PL,PR,LIMIT) 
          CALL INTERPOLATION(DD(I,J-1:J+2),DL,DR,LIMIT) 
          CALL INTERPOLATION(UU(I,J-1:J+2),VXL,VXR,LIMIT) 
          CALL INTERPOLATION(VV(I,J-1:J+2),VYL,VYR,LIMIT) 
          CALL INTERPOLATION(HH(I,J-1:J+2),HL,HR,LIMIT) 
          CALL INTERPOLATION(GA(I,J-1:J+2),GAMML,GAMMR,LIMIT) 
	  ELSE 
          PL=PP(I,J)
	    PR=PP(I,J+1)
          DL=DD(I,J)
	    DR=DD(I,J+1)
          VXL=UU(I,J)
	    VXR=UU(I,J+1)
          VYL=VV(I,J)
	    VYR=VV(I,J+1)
          HL=HH(I,J)
		HR=HH(I,J+1)
          GAMML=GA(I,J)
		GAMMR=GA(I,J+1)
	  END IF 
      END IF                                                      
      IF(DIRECTION==1) THEN
        IF(I==0) THEN
          COS=-AREA_W(2,1,J)
          SIN=-AREA_W(3,1,J)
        ELSE
          COS=AREA_E(2,I,J)
          SIN=AREA_E(3,I,J)
        ENDIF
      ELSE IF(DIRECTION==2) THEN
        IF(J==0) THEN
          COS=-AREA_S(2,I,1)
          SIN=-AREA_S(3,I,1)
        ELSE
          COS=AREA_N(2,I,J)
          SIN=AREA_N(3,I,J)
        END IF
      END IF
      CL=DSQRT(GAMML*PL/DL)  
      CR=DSQRT(GAMMR*PR/DR)
      CI=.5D0*(CL+CR)  ! SQRT(CL*CR)
      VNL=VXL*COS+VYL*SIN
      VNR=VXR*COS+VYR*SIN
      MAL=VNL/CI 
      MAR=VNR/CI 
      IF(DABS(MAL)<1.D0)  THEN
        B1=.25D0*(MAL+1.D0)**2+.125D0*(MAL**2-1.D0)**2
        G1=.25D0*(MAL+1.D0)**2*(2.D0-MAL)+3.D0/16*MAL*(MAL**2-1.D0)**2
      ELSE
        B1=.5D0*(MAL+DABS(MAL))
        G1=B1/MAL
      ENDIF
      IF(DABS(MAR)<1.D0)  THEN
        B2=-.25D0*(MAR-1.D0)**2-.125D0*(MAR**2-1.D0)**2
        G2=.25D0*(MAR-1.D0)**2*(2.D0+MAR)-3.D0/16*MAR*(MAR**2-1.D0)**2 
      ELSE
        B2=.5D0*(MAR-DABS(MAR))
        G2=B2/MAR
      END IF
      MAI=B1+B2
      VI=MAI*CI
      IF(VI>=0.D0) THEN
        FI(1)=VI*DL
        FI(2)=VI*DL*VXL+(G1*PL+G2*PR)*COS
        FI(3)=VI*DL*VYL+(G1*PL+G2*PR)*SIN
        FI(4)=VI*(DL*HL+.5*DL*(VXL**2+VYL**2))
      ELSE
        FI(1)=VI*DR
        FI(2)=VI*DR*VXR+(G1*PL+G2*PR)*COS
        FI(3)=VI*DR*VYR+(G1*PL+G2*PR)*SIN
        FI(4)=VI*(DR*HR+.5*DR*(VXR**2+VYR**2))
      END IF 
      RETURN
      END SUBROUTINE AUSM
!----------------------------------------------------------------------C 
! PURPOSE:                                                             C
!----------------------------------------------------------------------C
      SUBROUTINE INTERPOLATION(POINT,LL,RR,FLAG)
      IMPLICIT NONE 
      REAL*8:: POINT(4),A0,A1,A2,LL1,LL2,LL,RR,R1,R11,R2,R22,K0,LR1,LR2
      REAL*8,EXTERNAL::MINMOD,VAN_LEER,DOUBLE_MINMOD,SUPERBEE,                  VAN_ALBADA
      INTEGER::FLAG
      K0=1.D0/3.D0
      A1=POINT(3)-POINT(2)
      A0=POINT(2)-POINT(1)
      A2=POINT(4)-POINT(3)
      R1=A1/(A0+1.D-12)
      R2=A1/(A2+1.D-12)
      SELECT CASE(FLAG)
      CASE(1) 
        LL1=MINMOD(R1) 
	  LL2=MINMOD(1.D0/R1)
        LR1=MINMOD(R2) 
	  LR2=MINMOD(1.D0/R2)
      CASE(2)
        LL1=VAN_LEER(R1) 
	  LL2=VAN_LEER(R11)
        LR1=VAN_LEER(R2) 
	  LR2=VAN_LEER(R22)
      CASE(3)
        LL1=SUPERBEE(R1) 
	  LL2=SUPERBEE(R11)
        LR1=SUPERBEE(R2) 
	  LR2=SUPERBEE(R22)
      CASE(4)
        LL1=VAN_ALBADA(R1) 
	  LL2=VAN_ALBADA(R11)
        LR1=VAN_ALBADA(R2) 
	  LR2=VAN_ALBADA(R22)
      CASE(5)
        LL1=DOUBLE_MINMOD(R1) 
	  LL2=DOUBLE_MINMOD(R11)
        LR1=DOUBLE_MINMOD(R2) 
	  LR2=DOUBLE_MINMOD(R22)
      END SELECT
      LL=POINT(2)+0.25D0*((1.D0-K0)*LL1+(1.D0+K0)*LL2*R1)*A0 
      RR=POINT(3)-0.25D0*((1.D0-K0)*LR1+(1.D0+K0)*LR2*R2)*A2 
      RETURN
      END SUBROUTINE INTERPOLATION
!----------------------------------------------------------------------C 
! PURPOSE:                                                             C
!----------------------------------------------------------------------C
      FUNCTION VAN_LEER(R)
      REAL*8:: VAN_LEER,R
        VAN_LEER=(R+ABS(R))/(1.D0+ABS(R))
      RETURN
      END FUNCTION VAN_LEER
!----------------------------------------------------------------------C 
! PURPOSE:                                                             C
!----------------------------------------------------------------------C
      FUNCTION VAN_ALBADA(R)
      REAL*8:: VAN_ALBADA,R
      VAN_ALBADA=(R**2+R)/(1.D0+R**2)
      RETURN
      END FUNCTION VAN_ALBADA
!----------------------------------------------------------------------C 
! PURPOSE:                                                             C
!----------------------------------------------------------------------C
      FUNCTION MINMOD(R)
      REAL*8:: MINMOD,R
      IF(R>0.D0) THEN
        MINMOD=MIN(R,1.D0)
      ELSE
        MINMOD=0.D0
      END IF
      RETURN
      END FUNCTION MINMOD
!----------------------------------------------------------------------C 
! PURPOSE:                                                             C
!----------------------------------------------------------------------C
      FUNCTION DOUBLE_MINMOD(R)
      REAL*8:: DOUBLE_MINMOD,R
      IF(R>0.D0) THEN
        DOUBLE_MINMOD=MIN(2*R,1.D0,(1D0+R)/2)
      ELSE
        DOUBLE_MINMOD=0.D0
      ENDIF
      RETURN
      END FUNCTION DOUBLE_MINMOD
!----------------------------------------------------------------------C 
! PURPOSE:                                                             C
!----------------------------------------------------------------------C
      FUNCTION SUPERBEE(R)
      REAL*8:: SUPERBEE,R
      SUPERBEE=MAX(0.D0,MIN(2*R,1.D0),MIN(R,2.D0))
      RETURN
      END FUNCTION SUPERBEE